Draft version February 2, 2008 

Preprint typeset using 1^1^^ style cmulatcapj v. 10/10/03 



EVOLUTION OF MASSIVE BLACK HOLE BINARIES 
Junichiro Making ^ and Yoko Funato^ 

Draft version February 2, 2008 

ABSTRACT 

We present the results of large-scale iV-body simulations of the stellar-dynamical evolution of massive 
black hole binaries at the center of spherical galaxies. We focus on the dependence of the hardening rate 
on the relaxation timescale of the parent galaxy. A simple theoretical argument predicts that a binary 
black hole creates a "loss cone" around it. Once the stars in the loss cone are depleted, the hardening 
rate is determined by the rate at which field stars diffuse into the loss cone. Therefore the hardening 
timescale becomes proportional to the relaxation timescale. Recent iV-body simulations, however, have 
failed to confirm this theory and various explanations have been proposed. By performing simulations 
with sufhciently large iV (up to 10^) for sufficiently long time, we found that the hardening rate 
does indeed depend on N . Our result is consistent with the simple theoretical prediction that the 
hardening timescale is proportional to the relaxation timescale. This dependence implies that massive 
black hole binaries are unlikely to merge within a Hubble time through interaction with field stars 
and gravitational wave radiation alone. 

Subject headings: black hole physics — galaxies:interactions — galaxies: nuclei — methods: N-body 
simulations — stellar dynamics 



1. INTRODUCTION 

1.1. Massive black hole binaries and the central 
\ structure of ellipticals 

The possibility of the formation of massive black hole 
binarie s in the cores of e lliptical galaxies was first pointed 
out by 'BcEelman et al. (1980|, hereafter BBR). At that 
time, the possibility of such events was purely theoreti- 
cal, derived from two hypothesis. One is that QSOs are 
driven by massive central black holes, with masses of the 
order of IO^Mq or larger. The second is that most el- 
liptical galaxies, which might contain such massive black 
holes, are formed through merging of two galaxies. If 
both of the progenitor galaxies contain massive black 
holes, these black holes sink to the center of the merger 
product through dynamical friction and form a binary. 

Though this scenario was purely theoretical at the 
time first proposed, a considerable amount of observa- 
tional support has been provided in the last two decades. 
First of all, there is now plenty of evidence that many, if 
not m ost, giant ellipticals con tain massive central black 
holes l|Magorrian et al.lll998j) . Also, it has been sug- 
gested that the black hole mass Mbh shows a tight cor- 
relation with the spheroidal mass and the central veloc- 
ity di spersion ijGebhardt et al.l2000l : lFerrarese &: Merritti 
I2000D . The most straightforward and natural explana- 
tion of the observed correlations is that massive galax- 
ies are formed by merging of less massive galaxies, 
and that the central black als o grows through merging 
l|Kauffmann fc HaehneIill2nnflD . 

This merging scenario has an additional important ad- 
vantage: it nicely explains the observed structure of the 
central region of massive elliptical galaxies. In the 1970s 
and 1980s, observations revealed that large ellipticals had 
large "cores" , with a good linear correlation between the 
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size of the core and the size of the galaxy. At that time, 
this correlation was thought to be a very strong counter 
evidence against the merger hypothesis. The larger core 
size of the larger galaxy implies that the central phase 
space density is lower for larger galaxies. However, if 
larger galaxies had been formed by coUisionless merging 
of smaller galaxies, the central phase space density would 
have been roughly conserved, and therefore the core size 
would also have been roughly conserved. 

Strictly speaking, the phase space argument only guar- 
antees that the central phase space density would not 
increase, and does not imply that it is con served. Never- 
theless, numerical simulations of mer gers (iFarouki et alJ 
ri98^ White m§, lOkumuraet in 119911) all demou- 
strated that the core size would not increase significantly 
through merging. And if some dissipational process im- 
plicit in the dynamics of gas clouds is taken into account, 
most likely the central density would increase, rather 
than decrease. Thus, it seemed very difficult to explain 
the observed correlation between the core size and the 
size of the galaxy within the framework of the merger 
h ypothesis. 

lEbisuzaki et alJ l)1991|) proposed that the formation of 
a black hole binary might solve this difficulty. If both 
progenitors in a collision between two galaxies contain 
central massive black holes, these black holes would form 
a binary as suggested by BBR. The back reaction of the 
sinking of black holes through dynamical friction and the 
subsequent hardening of the black hole binary cause a 
heating of the stars in the core, leading to an expansion of 
the core. Using iV-body simulations of spherical galaxies 
with and without central black holes, they demonstrated 
that the core size, defined as the density-weighted dis - 
tance from the density center ijCasertano fc Hu'3ll985(l . 
increases in the case of a merger of two galaxies with 
central black holes. 

In the mid-1990s, high-resolution observations by HST 
revealed that the "cores" of the giant elliptical galaxies 
are not really cores with a flat volume density, but rather 
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very shallow cusps in which the density continues to rise 
toward the center as p oc r~", where p is the volume 
luminosity density and r is the distance from the center, 
with power index a = 0.5 ~ 1 (Gebhardt et al. 1996; 
iBvun et al. 19961. 

Around the same tim e, the completion of the GRAPE- 
4 l|Makino et al.llT997j) made it possible to perform A^- 
body simulations of merging galaxies with central black 
holes using far more p articles than had been prev iously 
possible. For example, Ma kino fc Ebisuz aki (1996) per- 
formed a simulation of repeated mergers of galaxies with 
central black holes, where the final product of one merger 
simulation is used as the progenitor for the next simula- 
tion. They found that the central structure of the merger 
remnant converges to a unique profile, in the form of a 
central cusp around the black hole with a slope of ap- 
proximately —0.5; the total mass of the stars in the cusp 
region is com parable to t he mass of the black hole binary. 
Furthermore. INakano fc Makino ( 1999a b) showed, by a 
combination of A^-body simulations and analytic argu- 
ments, that this shallow cusp can be explained by the 
fact that the distribution function of stars has a lower 
cutoff energy Eq. If there is a lower cutoff in energy dis- 
tribution, one can show in a simple derivation that the 
density profile close to the central black hole must have 
a slope with a power of —0.5. Such a cutoff is therefore 
naturally formed when black holes sink to the center and 
form a binary. 

In this theory, the radius of the cusp region relative 
to the effective radius of the galaxy is proportional to 
the mass of the black hole relative to that of the par- 
ent galaxy. Thus, the correlation between the spheroida l 
mass and the black hole mass l|Magorrian et al.l 11998^1 
nicely explains the correlation betw een the cusp radius 
and t he effective radius. Indeed, iMilosavlievic et all 
(|2002tl estimated "the mass ejected from the central 
cusp" of observed ellipticals, and found that it correlates 
well with the black hole mass. The correlation between 
the spheroidal mass and the black hole mass itself is thus 
explained by the merger scenario. 

1.2. The Fate of a black hole binary 

As we reviewed in the previous section, a merger of 
galaxies with central black holes can account for the ob- 
served characteristics of the central structure of elliptical 
galaxies quite well. However, we are still faced with one 
fundamental question: what will happen to the binary 
black hole? 

Once formed, a binary black hole further hardens 
through interactions with nearby stars, much in the same 
way as hard bi nari es in globul ar clusters continue to 
harden ISmV^ fl987t iHeggie fc l20nl . If we can 
assume that the stellar density around the black holes 
is roughly constant in time, the hardening rate of the 
binary, dEb/dt, is also constant. 

However, there are two significant differences. One is 
that the black hole binary is far more massive than the 
field stars. Therefore, it stays at the center of the galaxy 
with only a very small random velocity. The second is 
that a galaxy contains far more stars than globular clus- 
ters do, and therefore the central two-body relaxation 
time is much longer for a galaxy than for a globular clus- 
ter. BBR argued that because of these two differences, 
stars which can interact with the binary will sooner or 



later be depleted, and the hardening rate of the black 
hole binary will drop to practically zero. Once the nearby 
stars (stars with sufficiently small angular momentum to 
approach to the black hole binary) are scattered into dif- 
ferent orbits, in other words, after the "loss cone" is de- 
pleted, the hardening rate drops to a quite small value. 
In this stage, the hardening rate is determined by the 
timescale in which the loss cone is refilled by thermal re- 
laxation. Thus, in real ellipticals with very long central 
relaxation times, the timescale of refilling of the loss cone 
is very long and therefore the hardening timescale of the 
binary would be much longer than the Hubble time. 

Whether or not this loss-cone depletion actually occurs 
has been one of the key issues in the theoretical and nu- 
merical stu dies of m assive black hole binaries. E arly sim- 
ulations (jMakino e t al. 1993; Mikko la fc Valtonen 199^ 
could not cover a large enough range for the number 
of particles (and thus for the relaxation time) to see 
the corr e spond ing change in the hardening timescale. 
IMakinol l|1997D performed merger simulations with a 
number of particles A^ in the range of 2K to 256K, 
and found that the hardening timescale showed a clear 
increase, for increasing A^. From numerical results, 
Makino estimated that the hardening timescale was 
proportional to N^^^. This dependence is notably 
weaker than BBR's theoretical prediction that the hard- 
ening timescale sho uld be proportional t o the relax- 
ation time, or A^. IQuinlan fc HernauistI ljl997r per- 
formed similar calculations with up to 200K stars, and 
found a qualitatively similar dependence of the harden- 
ing timescale on the number of stars. However, they 
observed that the evolution of the black hole binary 
was almost the same for the runs with lOOK stars 
and the run with 200K stars, and concluded that for 
large enough A^ the hardening rate would become in- 
dependent of A^. Xlhatterjce et al. (2003) rep eated sim- 
ilar ca lculations as performed bv lOuinlan fc Hernauis^ 
ljl997fl . but using up to 400K stars, and reached the 
same conclusion Both lOuinlan fc Hernauisj l)1997|) and 
iChatteriee et al.' ('2003') used a combination of the SCF 
technique (Hernauist fc Ostriker 1992) and direct calcu- 
lation, where inte ractions betwee n field stars are treated 
usmg SCF, while IMakinol |T997(l relied on a fully direct 
calculation of all forces using the GRAPE-4. Though the 
SCF approximation docs not eliminate two-body relax- 
ation (Hernauist & Barnes 1990; Hernauist & Ostrike^ 
Il992j) . it may have affected the dependence of the re- 
laxat ion timescale on A^ in a com plex way. 

Mil osavlievic fc MerrittI ()2001f l performed simulations 
of mergers with central black holes, by again using a com- 
posite method, but here composite in the time domain 
and not in the spatial domain. Before the formation of 
the black hole binary, the y employed a pa rallel treecode 
with individual timesteps ijSpringel et al.l2 001'). Just be- 
fore the binary formation, they then switched over to a 
direct calculation using a general-purpose parallel com- 
puter. The use of a general-purpose parallel computer 
for direct calculations limited the number of particles to 
32K and less, but their simulation is unique in that they 
started from a self-consistent model with a central cusp 
of slope —2. They did not find any dependence of the 
hardening rate on the number of stars. 

To summarize, there is no single accepted view for 
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the evolution of massive black hole binaries in the cen- 
ter of galaxies. The simple analytic theory (BBR) pre- 
dicts that the hardening timescale should be determined 
by the timescale for refilling the loss-cone, and there- 
fore should be proportional to the relaxation time, or 
roughly speaking, the number of stars in the system. 
The re sults of numerical studies ran ge from no depen- 
dence llMilos avlievic fc MerrittI 120011) to a dependence 

DC Ari/3 ( Makii3[lM3)- 

In this paper, we give a clear and decisive answer to 
the question whether the hardening timescale depends 
on the number of particles used in the simulation, and 
if so, in what way. In section 2 we describe the initial 
model and the numerical method used. In section 3 we 
give the result. Section 4 contains the discussion. 

2. NUMERICAL METHOD AND INITIAL MODELS 
2.1. Numerical method 

We performed A^-body simulations of galaxies with 
massive black holes. For all calculations, we used a 
program with direct force calculation and fourth-order 
Hermite i ntegrator with individual (block) timestep 
(^akino fc Aarsethlll992(l . For gravitational interaction 
between field stars, we apply the usual Plummer soft- 
ening. The size of the softening is described in section 
12.21 The gravitational interactions between black holes 
and that between black holes and field particles are cal- 
culated with very small softening, ebh — 10~^. Since we 
do not apply any regularization technique, we need to 
guarantee that the gravitational force does not diverge. 

For calculation of gravitational forces from field parti- 
cles (to both the field particles and black h oles), we used 
the s pecial-purpose hardware GRAPE-6 ljMakino et al.l 
l2002tl . The calculation of the forces from black holes 
was done on the host comp uter to maintain suffic i ent ac - 
curacy, in a way similar to IMakino fc Ebisuzakil l|1996fl . 
The relative accuracy of the pairwise force calculated 
with "High accuracy" versions of GRAPE hardwares 
(GRAPE-2, 4 and 6) is roughly single precision (24-bit 
mantissa). This accuracy is usually sufficient, since the 
errors of the forces from many particles partially cancel 
each other. However, single-precision is not quite enough 
for forces from black hole particles. The reason is that 
some particles, and most notably the black hole particles 
themselves after they formed a binary, orbit around a 
black hole particle a number of times, essentially feeling 
only the force from one black hole particle. In this case, 
the round-off error accumulates, and the error in the to- 
tal energy becomes alarmingly large. Moreover, all the 
errors are generated from the black hole binary and stars 
which are close to them, the behavior of which we are in- 
terested in. In order to guarantee sufficient accuracy, 
we chose to calculate all gravitational forces from black 
hole particles on the host computer, using full double- 
precision (53-bit mantissa) numbers. 

No relativistic effect was taken into account. Our cal- 
culation is purely Newtonian. We do not model the ac- 
cretion to black holes or collisions between stars, since 
their cross sections are small. 

2.2. Initial Models 

In this paper, we consider a simple model, where we 
place two massive point-mass particles in a spherical 



galaxy. For our standard set of runs, the initial galaxy 
model is a King model with nondimensional central po- 
tential Wq — 7. We use Heggie units, where the mass M 
and the virial radius Ry of the initial galaxy model and 
gravitational constant G are all unity. In these units, the 
binding energy of the initial galaxy is E = —1/4. 

The mass of the black hole particles is Mbh = 0.01. 
They are initially placed at (±0.5,0,0) with velocity 
(0, ±0.1,0). Thus, they are initially outside the core, 
at the apocenter of nearly radial orbits. 

We varied the number of particles N from 2,000 to 
1,000,000. To investigate the possible dependence of the 
results on the softening length, we tried three different 
choices for the softening length: (a) e — 0.01, (b) e = 
0.01/(iV/2000)i/3, (c) e = 20/A^. AU give the same e = 
0.01 for N ~ 2000, but they have a different dependence 
on N. 

The largest calculation (1 million particles for up to 
t = 300) took about one month on a single-host, 4- 
processor-board GRAPE-6 system with a peak speed 
of 4 Tflops. The total number of individual timesteps 
was 1.2 X 10^^. In other words, the (harmonic) average 
timestep is around 2.5 x 10~^. In comparison, the total 
number of block timesteps is 2.1 x 10*. Thus, the typical 
timestep size of the particle with the smallest timestep 
(generally the black hole particles themselves or particles 
close to them) is 1.4x 10^^, more than 1000 times smaller 
than the average stepsize. Without the use of individual 
timesteps, this calculation would have required a fraction 
of a Petaflops-year, rather than 1/3 Teraflops-year. 

For all calculations, the total energy is conserved to 
better than 0.1%, or better than 1% of the binding energy 
of the BH binary. We did several test calculations with 
both higher and lower accuracy criteria, but found no 
systematic difference in the final results. 

3. RESULTS 

3.1. The N -dependence 

Figure n shows the time evolution of the specific bind- 
ing energy (per unit of reduced mass) of the black hole 
binary. The relation between the semi-major axis a and 
the binding energy is 



The relation between the orbital velocity of the black 
holes and Ejj is given simply by 

VBH = Vc^/\E'b\, (2) 

where Vc is the three-dimensional velocity dispersion of 
the field stars, in a central region of the galaxy, chosen 
to be large enough not to be affected significantly by the 
black holes. Thus, if we interpret the host galaxy as a 
galaxy with a velocity dispersion of 300 km/s, the black 
hole binary with Eb = —1 has an orbital velocity of 300 
km/s. 

From figure ^ it is clear that the evolution timescale 
continues to depend on N for the entire range of N for 
which we performed our simulations, with no indication 
whatsoever of even an onset toward convergence. 

Figure El shows the early evolution. Initially the bind- 
ing energy shows large oscillations, simply because the 
black hole particles are not bound to each other but or- 
bit within the parent galaxy in highly eccentric orbits. 
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Fig. 1. — The evolution of the specific binding energy i?;, of the 
binary black hole. Curves give the results for A'^ = lOK, 20K, 50K, 
lOOK, 200K, 500K and IM (left to right). For all calculations, the 
softening length is e = 0.01. 




Fig. 2. — Same as figure^ but only the early evolution is shown. 
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Fig. 3. — The hardening rate /3 = |Ai?b/At| plotted as a function 
of the number of particles A'^. Dotted, dash-dotted, dashed and 
solid curves denote the values measured from the time at which 
\E}j\ reached 1, 3, 5 and 7, respectively, until the time at which 
has increased by an additional amount of +0.5. 
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Fig. 4. — The slope 7 of the dependence of the hardening rate /3 
on the number of particles A'^, plotted as a function of the binding 
energy -Ej o '^^ which /3 is measured. 



As their orbits shrink through dynamical friction, the 
ampUtude of the osciUations becomes smaller, and even- 
tually the two black holes become bound. We can see 
that the early evolution of the black hole binary, before 
the specific binding energy reaches —1, is almost indepen- 
dent of N . However, after E}, reaches —2, the evolution 
timescale shows a strong dependence on the number of 
particles. 

To quantitatively evaluate the dependence of the hard- 
ening rate on the number of particles, we calculate the 
hardening rate (3, defined as 

(3=\^Eb/M\. (3) 

Here, At — ti—tQ where a-nd ti are the times at which 
Eb reached the values E^^o and E^fl + AE^, respectively. 
We use AEii = — 0.5 for all values of E^ q. FigureOshows 
the result, for Eb^ = -1,-3,-5,-7. When Eb = — 1, 



the hardening rate is almost independent of N. How- 
ever, as the binary becomes harder, f3 decreases, and the 
decrease is larger for larger N. Thus, the hardening rate 
/? for large values of \Eb\ shows a strong dependence on 
the number of particles N. 

This result is exactly what is expected from the sim- 
ple loss-cone argument: after the loss cone is depleted, f3 
should be inversely proportional to the relaxation time, 
which is proportional to N. Note that we used a con- 
stant softening, for which the Coulomb logarithm does 
not depend on N. 

If the loss-cone argument would be correct, writing (3 oc 
would let 7 approach unity for large enough N. 
Figure 0] shows 7, calculated using the value of (3 for 
N = 20K and IM. The choice of these values of N is 
somewhat arbitrary, but as one can see from figure |31 
we obtained similar figures with different choices of N. 
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We could use a least-square fit, but decided against it 
since there is no obvious reason to assign equal weights 
to results with different N. 

We can see that 7 indeed increases as \Eb\ increases. 
Even at Ei, = —8, 7 has not converged to a final value. 
If we could extend the calculation to higher values of Eb, 
we would be able to determine whether or not 7 really 
approaches unity. Unfortunately, it would be too time 
consuming, even on a GRAPE-6, to further extend the 
calculations, given that the hardening rate is so small for 
large N. 

From the current simulations, we can safely conclude 
that the hardening rate /3 depends on the number of 
particles TV, and the power index of the dependence 7 
numerically obtained is larger than 0.7. The numerical 
result is consistent with the simple loss-cone argument, 
which predicts 7 = 1. 

3.2. Dependence on the softening length 

The results given in the previous subsection demon- 
strate clearly that the hardening rate depends on the 
number of particles. In this and following subsections, 
we will check whether this result is really reliable. First 
we look at the effect of the softening. The relatively large 
and constant softening used in our standard runs has the 
effect of suppressing two-body relaxation. In particular, 
in the core of the galaxy, the effective Coulomb logarithm 
might be very small, resulting in unphysical suppression 
of the relaxation effect. Since the timescale of the loss- 
cone refill is related to the relaxation timescale, it is cru- 
cial to express the relaxation with a reasonable accuracy. 
To test the effect of the softening, we performed sev- 
eral runs with a much smaller softening length. Figure |S1 
shows the result. For N = 2x 10^, reducing e by a factor 
of 100 resulted in only a small increase in the hardening 
rate. Furthermore, in the case of runs with N = 10^, the 
change in the softening has practically no effect on the 
hardening rate. Thus, we can safely conclude that the 
effect of the softening is, if any, sufhciently small that it 
does not affect the results in the previous section. 

3.3. Dependence on the initial models 

Before the loss cone is depleted, the hardening rate is 
expected to be proportional to the central density of the 
parent galaxy. After the loss cone is depleted, the hard- 
ening rate is determined by the timescale at which the 
stars close to the loss cone diffuse into the cone. Thus, 
here again, the hardening timescale is expected to de- 
pend on the initial central density. 

Figure ini shows the results of runs with different initial 
galaxy models. As expected, the hardening is faster for 
models with deeper central potential (higher central den- 
sity). However, for all runs the hardening rate depends 
on the number of particles. 

Figure [3 shows the growth rate /3 as a function of the 
initial central density of the parent galaxy, po^. In the 
early phase (Ef, — —1), /3 is roughly proportional to p for 
p < 10^. For higher p or for later phases the dependence 
is somewhat weaker, presumably because the black hole 
binary already has ejected nearby stars, thereby reducing 
the central density. 

3.4. Loss cone depletion and wandering of the binary 
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Fig. 5. — The evolution of specific binding energy . The curves 
denote, from top to bottom, AT = 10^, e = 0.01 (sohd), e = 0.00045 
(dashed), AT = 2 X 10^,e = 0.01 (soHd), t = 0.0001 (dashed), and 
e = 0.00001 (dotted). 
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Fig. 6. — The evolution of specific binding energy Ei, for runs 
with different initial galaxy models. Solid, dashed and dotted 
curves display the results for Af = lO'', 2 X 10^^, and 2 X 10*. The 
left three curves are for a King model with ipo = H, while the right 
three curves are for ipo = 7. 



In the previous subsections we have seen that the be- 
havior of the hardening rate is consistent with the loss 
cone argument. In this subsection, we directly investi- 
gate whether or not the loss cone is actually depleted. 

Figure |H1 shows the distribution of particles in the 
{E, J) plane, where E is the specific energy and J is the 
specific total angular momentum. We use the coordinate 
origin as the reference point for the angular momentum. 
We also tried to use the center of mass of the black hole 
binary, but that resulted in practically indistinguishable 
figures. Here, we can clearly see that the number of par- 
ticles with J < 0.01 is more and more depleted as time 
proceeds. At T = 10, only particles with low energy 
{E < —2) are affected. However, the depletion reaches 
higher energy levels as time proceeds, and by T = 80 
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Fig. 7. — The growth rate j3 = \dE},/dt\ plotted as a function 
of the initial central density po,0 of the parent galaxy. Different 
curves denote the growth rate at different values of | | . 
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Fig. 8. — The distribution of particles in the {E,J) plane at 
times T = 10, 20, 40 and 80 (from top left to bottom right). The 
number of particles is 10® . 



only stars with nearly zero energy are left with low an- 
gular momentum. Note that there were no particles with 
nearly zero energy at T = 10. These barely bound parti- 
cles have been kicked to high-energy, long-period orbits 
by the binary black hole. 

Figure El shows the radial density profiles for the same 
snapshots as used for Figure |H1 Here we took the center 
of mass of the black hole binary as the center of the 
coordinate. Though the loss cone depletion is clearly 
visible in the distribution of particles in the {E, J) plane, 
the density profiles show no clear sign of the existence of 
a loss cone. This result is not at all surprising. If the 
distribution of particles in the {E, J) plane is the same 
as that of the initial King model, adding the central black 
hole potential would result in a density cusp with a slope 
of -0.5 pakano k Makino 1999b). Thus, the fact that 
the density docs not increase toward the center actually 
means that the particles close to the black hole binary 
are strongly depleted. 
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Fig. 9. — Radial density profiles at times 10, 20, 40 and 80 for 
the run with N = 10^ . 
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Fig. 10. — The distribution of particles in the {E,J) plane for 
runs with 10* (left) and 10^ (right) particles, at t = 60. 



Figure ^1 shows the distribution of particles in the 
{E, J) plane for runs with smaller number of particles. 
The distribution for = lO"' does not show any sign of 
depletion, while the distribution for N — 10^ shows a 
weak indication of depletion. It also show an enhance- 
ment of nearly zero-energy stars with small angular mo- 
mentum. 

Note that figure |H1 shows that the periastron distance 
of the particles that are depleted, a measure for the ef- 
fective radius of the loss cone, is 0.01 or larger. This 
radius is much bigger than the semi-major axis of the bi- 
nary, which is around 0.001 at T = 80. The reason why 
the effective size of the loss cone is much bigger than the 
semi-major axis of the black hole binary is the wandering 
of the binary. Figure lTTl shows the time variation of the x 
coordinate of the center of mass of the binary. The typi- 
cal wandering distance is roughly proportional to N~^^'^, 
as theoretically p redicted and demonstrated by previous 
numerical work ijMakinol Il997t iMilosavlievic &: Merritd 
120011) . For N = 10^, the typical wandering distance is 
around 0.01, roughly consistent with the size of the loss 
cone. 

4. DISCUSSION 

4.1. Comparison with previous works 

As we summarized in the introduction, several re- 
searchers have performed iV-body simulations of the evo- 
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Fig. 11. — The time variation of the x component of the center of 
mass of the black hole binary. The number of particles is 10* (top), 
10^ (middle), and 10^ (bottom). Note that the vertical axis is scaled 
in proportion to l/\/~N. 



lution of massive black hole binaries in the center of a 
galaxy, and nobody has obtained a result which could 
be interpreted as being consistent with the simple loss- 
cone argument bv lBeeelman et all l|198ClD . In the present 
work, however, we have obtained a result that is consis- 
tent with the loss cone argument. In order to under- 
stand the cause of the discrepancy, we first discuss our 
own previous work (iMak ina.,1997 |). and then work by oth- 
ers lIQuinlan fc Hernauistill997t iMilosavlievic fc Merritti 
12001^. 'Chatteriee et al.ll2003() . 

[Makino ( 1997) obtained a value for the slope 7 of the 
dependence of the binary hardening rate P on the number 
of particles iV of 7 '--^ 1/3, from merger simulations with 
N ~ 2, 048 to 262,144. If we compare his figure 1 and our 
figurcs^andl^ the reason why the value of 7 was small is 
obvious: the value of the binding energy at which /? was 
measured was simply too small. Therefore f3 had not yet 
reached the value determined by the relaxation timescale. 
As shown in figures |3I and ^ 7 increases as \Eb\ increases. 
In hindsight, it is clear that previous simulations did not 
c over a long enough tim e . 

iQuin lan fc HernauistI l)1997j) performed simulations 
very similar to those presented here. They found that 
/3 depended on N, for N < 10^. However, the hard- 
ening rate was practically the same for the runs with 
N = 10^ and iV = 2 x 10^. They explained this result 
as follows. In their calculation, they had used the SCF 
method to evaluate the gravitational interaction between 
field stars. Therefore the two-body relaxation had been 
strongly suppressed, and the dependence of the harden- 
ing rate on the number of particles should come only 
through the wandering (Brownian motion) of the black 
hole binary. The random velocity of the black hole bi- 
nary would be proportional to V^iV. However, the dis- 
tance the wandering covers would not become arbitrary 
small, since the black hole binary would deplete the "loss 
cone" and create a kind of vacuum around it, the black 
hole binary would not feel a restoring force as long as it 
remained in the vacuum. 



The a bove explanation loo k s pla usible, but unfor- 
tunately IHuInla^^Simanial ^^1997^ gave no evidence 
that the distance covered by the black hole binary be- 
came independent of N. Our simulations demonstrated 
that the excursion distance is proportional to 1/ ^/N, for 
up to 10^ particles. 

An additional complication with their calculations is 
that it is difhcult to estimate the effects of two-body 
relaxation and wandering, because they used a mass for 
the particles which depended on their initial angular mo- 
menta (or, precisely, the periastron distances). They 
used this trick to increase the mass resolution near the 
center of the galaxy. However, because of this trick, the 
typical mass of stars that interact with the black hole bi- 
nary would become larger as the black hole binary kicks 
out more and more of the nearby (low-mass) particles. 
This implies that the strength of the Brownian motion 
depends on how much mass has been kicked out of the 
core. Also, even though the effect of direct two-body en- 
counters of field particles is suppressed, these particles 
can still indirectly exchange energy and angular momen- 
tum since all field stars directly interact with the black 
hole binary which has a finite random velocity. Even 
distant massive particles do interact significantly with 
the black holes. Thus, the way the two-body relaxation 
scales is quite difficult to evaluate for their scheme. 

Quinlan and Hernquist tried one simulation in which 
the center of mass of the binary is fixed at the origin 
of the coordinates, and found that the hardening rate 
dropped dramatically. They argued that this drop is ev- 
idence for wandering being important. However, there 
is another, equally plausible explanation. In their SCF 
expansion, they incorporated only spherically symmetric 
terms. Thus, when they fixed the center of mass of the 
binary to the origin, the gravitational potential calcu- 
lated became strictly spherically symmetric, except for 
the multipole moment of the binary black hole. Thus, 
the angular momentum of field stars with periastron dis- 
tance larger than the semi-major axis of the binary is 
conserved. In other words, the loss cone, once depleted, 
can never be refilled. Once they allow the center-of-mass 
motion of the binary to occur, however, the gravitational 
potential is no longer spherically symmetric, and angu- 
lar momenta (as well as energies) of field particles change 
through interactions with the center-of-mass motions of 
the binary. This process effectively works as a relaxation 
mechanism and causes the refilling of the loss cone. 

To summarize, though lOuinlan fc Hernauisti l|1997fl 
observed that in their calculation the hardening rate 
became independent of N for large N, it is some- 
what difhcult to generalize their result, obtained with 
a spherically-symmetric potential expansion code and 
radius-dependent particle mass, to real A^-body system 
o r real elliptical ga l axies. 

iChatteriee et al.l |2003') performed calculations similar 
to that presented in Quinlan & Hernauist (1997), using 
essentially the same calculation code. However, they did 
use a few non-spherical terms in their potential expan- 
sion, and they used equal-mass particles instead of par- 
ticles having radius-dependent mass. Thus, their results 
can be directly compared with our results without much 
complication. In the text, they stated that the hard- 
ening rate became constant for the value of N around 
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2 ~ 4 X 10^ For their 400K run, Mbh = 0.00125, which 
is 1 /8 of the value we used. Unfortunately, the simulation 
was stopped before we would expect to to see whether 
the loss-cone depletion would or would not occur. 

Unlike the p reviou s two papers, 

iMilosavlievic & Merritti l|200lD calculated the actual 
me rger of t wo g alaxies with central black holes, 
as 'Maki nol lll997l) di d. The difference between 
IMilosavlie vic & MerritH l|20fllD and iMakinol l|1997f) 
is that the former started with a galaxy model with 
a density cusp with p oc around the central black 
hole, while the latter used a King model with a finite 
core as the initial model. For simulations of mergers 
of elliptical galaxies with massive central black holes 
(IO^Mq or larger), a galaxy model with finite-size 
core would be more appropriate, since the shallow 
"cusp" of such large ellipticals corresponds to a cutoff 
of the distribution function at finite energy. The stellar 
distribution in the central region of spiral galaxies 
is consistent with a cusp of p cx r~^. Therefore, for 
simulations of mergers of spiral galax i es, th e initial 
model used by IMilosavlievic fc MerrittI l)2001|) may be 
more appropri ate. 

Milosavlicvi c fc Merritt) l)2001|) performed three runs 
with 8K, 16K and 32K stars, and found that the harden- 
ing rate was independent of the number of particles. As 
discussed in their paper, this result was stemmed from 
the fact that the loss cone was not depleted in their sim- 
ulations, because of two reasons. The first is that the 
initial central density of their model was high, since they 
tailored the distribution function so that the progeni- 
tor galaxies had central cusps around the black holes. 
This is a situation very different from that used in other 
studies, where black holes were placed off-center in a sin- 
gle galaxy or placed at the centers of galaxies with rela- 
tively large cores. Since the initial stellar density around 
the black holes is very high, the binary initially hard- 
ens very rapidly. A binary with a small semi-major axis 
has a small interaction cross section, which implies that 
it takes a long time to deplete the nearby stars. The 
second reason is the relatively small number of particles 
employed in their simulation, which resulted in rather 
large random velocities. Therefore, Brownian motion of 
the binary covered a fairly large radius. Even in their 
largest calculation with 32K particles, the total mass of 
the stars which can enter the region covered by the Brow- 
nian motion of the black hole binary is much larger than 
the mass ejected by the binary. Thus, the very high 
initial central density and the relatively small number of 
particles conspired together to prevent the loss cone from 
bein g depleted. 

If iMilosavlievic fc Merritti l)200H) could have used a 
much larger number of particles, the wandering dis- 
tance would have shrunk, and the loss cone would have 
been formed early on. According to their own estimate 
(IMilosavlievic fc Mcrritt 200^), for iV > 2 x lO""^ the loss 
cone depletion becomes important. 

To summarize, there is no real dis crepancy among the 
results of full iV-body simulations. IMakinol l)1997j) ob- 
served weaker dependence of the hardening rate than 
obtained in the present study, simply because his cal- 
culations were not long enough. Milosavlicvi c fc MerrittI 
((2001,) found no dependence on A^, essentially because 
the number of particles they used was too small for the 



loss cone to be depleted. 

4.2. Merging timescale of massive black hole binaries in 

ellipticals 

As first suggested by BBR and confirmed by a number 
of foUowup works, if the hardening timescale of the black 
hole binary is proportional to the relaxation timescale 
of the parent galaxy, the evolution timescale of typical 
black hole binaries in elliptical galaxies exceeds the Hub- 
ble time by many orders of magnitude. In other words, 
binaries are unlikely to merge through encounters with 
field stars and gravitational wave radiation. 

Our results strongly suggest that the hardening 
timescale is indeed determined by the relaxation 
timescale, for large enough N and after the binary be- 
comes sufficiently hard. This implies that gravitational 
interaction with field stars is insufficient to let the binary 
merge. 

There are a number of alternative mechanisms that 
may lead to the merger of the black hole binary. If there 
is a significant amount of gas left at the center, or if gas is 
supplied from the disk during the merger event, it would 
certainly change the entire picture. However, in the case 
of a merger of two ellipticals, there is not much cold gas 
left in the resulting galaxy. In this case, the most likely 
outcome is that the binary, stuck at a certain semi-major 
axis, stays at the center of the loss cone. 

If the eccentricity of the binary goes up, the timescale 
of orbital evolution by gravitational wave radiation is 
reduced significantly. Roughly speaking, if the eccen- 
tricity reaches 0.9, a fair fraction of the binary black 
holes would merge in a time less than the Hubble time. 
Some of the early A'^-body simulations and scattering ex- 
periments have focused on this p ossibility l|Makino et alJ 
119931 iMikkola fc Valto'nenlll992() . However, the general 
consensus seems to be that the eccentricity does not 
change much during the hardening. 

This situation can change if there are more than two 
massive black holes. If the binary black hole has a long 
lifetime, it is quite natural to assume that some of the 
galaxies which contain binary black holes will undergo a 
further merger with another galaxy with a central mas- 
sive black hole or a binary. If we regard the black holes 
as point-mass particles interacting through Newtonian 
gravity, then with three (or more) black holes we ex- 
pect at most one black hole binary to be left in the 
galaxy, having ejected all othe r black holes by th e gravi- 
tational slingshot mechanism l)Saslaw et al .11 19 73) . How- 
ever, here the eccentricity effects might play an impor- 
tant role. Simple estimates assuming a thermal distri- 
bution of eccentricities ( Makino fc Ebisuzaki 1994) sug- 
gest that, during repeated three body interactions, the 
eccentricity of the binary can reach a very high value, 
resulting in quick merging through gravitational wave 
radiation. In principle, it is possible that multiple black 
holes form a stable hierarchical system, where evolution 
of the outer binary is halted because of the loss cone de- 
pletion. In this case, however, the inner binary would 
typically have a semi- major axis of orde 1/10 of that of 
the outer binary, and the gravitational wave radiation 
timescale of the inner binary would generally be short. 
Also, the Kozai mechanism could play a role in increasing 
the eccentricity of the inner binary (,Blaes et al 200211 . 

So far, our primary attention has been directed to su- 
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permassive black holes in large ellipticals, since observa- 
tional evidence for massive dark objects is strongest for 
large ellipticals. If LSB galaxies had central black holes 
(which would be more like intermediate-mass black holes 
than massive black holes), some of them must have expe- 
rienced major merger events and therefore are likely to 
have contained multiple black holes at some point. Here, 
again, the crucial question is whether or not these black 
holes can merge. Since the black hole masses are smaller, 
if the central region of LSB galaxies is dominated by nor- 
mal stars, loss-cone depletion is less effective, and black 
holes can become very tight. Indeed, if we regard the 
field particles in our simulations as normal stars with 
typical masses around a solar mass, we have performed 
simulations of BH binaries with masses up to 10^ M©, 
for which we have seen that the hardening rate is still 
pretty high. However, if the central region is dominated 
by CDM, as suggested by both theory and observations, 
the thermal relaxation time would be practically infinite 



and loss-cone depletion would occur instantly. In the case 
of LSB galaxies, it is unlikely that two black holes merge 
during triple interactions, because the central potential 
is shallow. All black holes would be ejected from the 
parent galaxy before they would get a chance to merge. 
Thus, unlike large ellipticals, LSB galaxies are unlikely 
to display something like a Magorrian relation or M^h-a 
relation. 
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